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Abstract 

The use of a finite mixture of normal distributions in model-based clustering 
allows to capture non-Gaussian data clusters. However, identifying the clusters from 
the normal components is challenging and in general either achieved by imposing 
constraints on the model or by using post-processing procedures. 

Within the Bayesian framework we propose a different approach based on sparse 
finite mixtures to achieve identifiability. We specify a hierarchical prior where the 
hyperparameters are carefully selected such that they are reflective of the cluster 
structure aimed at. In addition, this prior allows to estimate the model using stan¬ 
dard MCMC sampling methods. In combination with a post-processing approach 
which resolves the label switching issue and results in an identified model, our ap¬ 
proach allows to simultaneously (1) determine the number of clusters, (2) flexibly 
approximate the cluster distributions in a semi-parametric way using finite mixtures 
of normals and (3) identify cluster-specific parameters and classify observations. The 
proposed approach is illustrated in two simulation studies and on benchmark data 
sets. 
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1 Introduction 


In many areas of applied statistics like economics, finance or public health it is often 
desirable to find groups of similar objects in a data set through the use of clustering 
techniques. A flexible approach to clustering data is based on mixture models, whereby 
the data in each mixture component are assumed to follow a parametric distribution with 
component-specific parameters varying over the components. This so-called model-based 
clustering approach (Fraley and Raftery, 2002) is based on the notion that the component 
densities can be regarded as the “prototype shape of clusters to look for” (Hennig, 2010) 
and each mixture component may be interpreted as a distinct data cluster. 

Most commonly, a finite mixture model with Gaussian component densities is fitted to 
the data to identify homogeneous data clusters within a heterogeneous population. How¬ 
ever, assuming such a simple parametric form for the component densities implies a strong 
assumption about the shape of the clusters and may lead to overfitting the number of clus¬ 
ters as well as a poor classification, if not supported by the data. Hence, a major limitation 
of Gaussian mixtures in the context of model-based clustering results from the presence of 
non-Gaussian data clusters, as typically encountered in practical applications. 

Recent research demonstrates the usefulness of mixtures of parametric non-Gaussian 
component densities such as the skew normal or skew-f distribution to capture non-Gaussian 
data clusters, see Friihwirth-Schnatter and Pyne (2010), Lee and McLachlan (2014) and 
Vrbik and McNicholas (2014), among others. However, as stated in Li (2005), for many 
applications it is difficult to decide which parametric distribution is appropriate to char¬ 
acterize a data cluster, especially in higher dimensions. In addition, the shape of the 
cluster densities can be of a form which is not easily captured by a parametric distribution. 
To better accommodate such data, recent advances in model-based clustering focused on 
designing mixture models with more flexible, not necessarily parametric cluster densities. 

A rather appealing approach, known as mixture of mixtures, models the non-Gaussian 
cluster distributions themselves by Gaussian mixtures, exploiting the ability of normal 
mixtures to accurately approximate a wide class of probability distributions. Compared 
to a mixture with Gaussian components, mixture of mixtures models impose a two-level 
hierarchical structure which is particularly appealing in a clustering context. On the higher 
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level, Gaussian components are grouped together to form non-Gaussian cluster distributions 
which are used for clustering the data. The individual Gaussian component densities 
appearing on the lower level of the model influence the clustering procedure only indirectly 
by accommodating possibly non-Gaussian, but otherwise homogeneous cluster distributions 
in a semi-parametric way. This powerful and very flexible approach has been employed in 
various ways, both within the framework of finite and infinite mixtures. 

Statistical inference for finite mixtures is generally not easy due to problems such as 
label switching, spurious modes and unboundedness of the mixture likelihood (see e.g. 
Friihwirth-Schnatter, 2006, Chapter 2), but estimation of a mixture of mixtures model is 
particularly challenging due to additional identihability issues. Since exchanging subcom¬ 
ponents between clusters on the lower level leads to different cluster distributions, while the 
density of the higher level mixture distribution remains the same, a mixture of mixtures 
model is not identifiable from the mixture likelihood in the absence of additional infor¬ 
mation. For example, strong identihability constraints on the locations and the covariance 
matrices of the Gaussian components were imposed by Bartolucci (2005) for univariate data 
and by Di Zio et al. (2007) for multivariate data to estimate finite mixtures of Gaussian 
mixtures. 

A different strand of literature pursues the idea of creating meaningful clusters after 
having fitted a standard Gaussian mixture model to the data. The clusters are determined 
by successively merging components according to some criterion, e.g. the closeness of the 
means (Li, 2005), the modality of the obtained mixture density (Chan et ah, 2008; Hennig, 
2010), the degree of overlapping measured by misclassihcation probabilities (Melnykov, 
2016) or the entropy of the resulting partition (Baudry et ah, 2010). However, such two- 
step approaches might miss the general cluster structure, see Appendix E for an example. 

In the present paper, we identify the mixture of mixtures model within a Bayesian 
framework through a hierarchical prior construction and propose a method to simultane¬ 
ously select a suitable number of clusters. In our approach both the identification of the 
model and the estimation of the number of clusters is achieved by employing a selectively 
informative prior parameter setting on the model parameters. 

Our choice of prior parameters is driven by assumptions on the cluster shapes assumed 
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to be present in the data, thus being in line with Hennig (2010) who emphasizes that, 
“it rather has to be decided by the statistician under which conditions different Gaussian 
mixture components should be regarded as a common cluster”. This prior specification 
introduces dependence among the subcomponent densities within each cluster, by pulling 
the subcomponent means on the lower level toward the cluster center, making the cluster 
distributions themselves dense and connected. On the higher level, the prior is based on the 
notion that the cluster centers are quite distinct from each other compared to the spread 
of the clusters. The choice of the hyperparameters of this hierarchical prior turns out to be 
crucial in achieving identification and is guided by a variance decomposition of the data. 

Regarding the estimation of the number of clusters, a sparse hierarchical mixture of 
mixtures model is derived as an extension of the sparse finite mixture model introduced in 
Malsiner-Walli et al. (2016). There, based on theoretical results derived by Rousseau and 
Mengersen (2011), an overfitting Gaussian mixture with K components is specified where 
a sparse prior on the mixture weights has the effect of assigning the observations to fewer 
than K components. Thus, the number of clusters can be estimated by the most frequent 
number of non-empty components encountered during Markov chain Monte Carlo (MCMC) 
sampling. In this paper, rather than using a single multivariate Gaussian distribution, we 
model the component densities in a semi-parametric way through a Gaussian mixture 
distribution, and again use a sparse prior on the cluster weights to automatically select a 
suitable number of clusters on the upper level. 

Specifying a sparse prior on the weights is closely related to Bayesian nonparametric 
(BNP) Gaussian mixture models such as Dirichlet process mixtures (DPMs; Ferguson, 
1983; Escobar and West, 1995). The sparse prior on the cluster weights induces clustering 
of the observations, similar as for DPMs which have been applied in a clustering context 
by Quintana and Iglesias (2003), Medvedovic et al. (2004) and Dahl (2006), among others. 
The hierarchical mixture of mixtures model we introduce is similar to hierarchical BNP 
approaches such as the hierarchical DPM (Teh et ah, 2006). Very closely related BNP 
approaches are infinite mixtures of infinite Gaussian densities such as the nested DPM 
(Rodriguez et ah, 2008), the infinite mixture of infinite Gaussian mixtures (Yerebakan 
et ah, 2014), and species mixture models (Argiento et ah, 2014) which directly work on the 
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partition of the data. We discuss in Sections 2.4 and 3.1 similarities as well as differences 
between our approach and BNP models. 

We finally note that the implementation effort to estimate our model is moderate and 
standard MCMC methods based on data augmentation and Gibbs sampling (see Fruhwirth- 
Schnatter, 2006) can be used. Several approaches proposed in the literature can be used to 
post-process the MCMC draws in order to obtain a clustering of the data and also to allow 
for cluster-specific inference. For our simulation studies and applications we adapt and ex¬ 
tend the method suggested by Friihwirth-Schnatter (2006, 2011) which determines a unique 
labeling for the MCMC draws by clustering the draws in the point process representation. 

The rest of the article is organized as follows. Section 2 describes the proposed strat¬ 
egy, including detailed prior specifications, and relates our method to the two-layer BNP 
approaches in Rodriguez et al. (2008) and Yerebakan et al. (2014). Clustering and model 
estimation issues are discussed in Section 3. The performance of the proposed strategy is 
evaluated in Section 4 for various benchmark data sets. Section 5 concludes. 


2 Sparse hierarchical mixture of mixtures model 

2.1 Model definition 

Following previous work on hierarchical mixtures of mixtures, we assume that N obser¬ 
vations yi, i — 1,..., N of dimension dirn(y,;) = r are drawn independently from a finite 
mixture distribution with K components, 

K 

p{yi\®,v) = ^2vkPk(yi\0k), ® = (o 1 ,...,0 K ), (i) 

k= 1 

with each component distribution Pk{yi\@k ) being a mixture of L normal subcomponents: 

L 

Pk{yi\0k) = ^WkifuiyilHw s ki )• (2) 

i=i 

In order to distinguish the component distributions on the upper level from the Gaussian 
components on the lower level, we will refer to the former ones as “cluster distributions”. 
For clustering the observations based on Bayes’ rule, the cluster weights r] = (rft ,... ,r] K ) 
and the cluster densities Pk{yi\0k ) on the upper level (1) are relevant. 
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Since the number of data clusters is unknown and needs to be inferred from the data, 
we assume that (1) is an overfitting mixture, i.e. the specified number of clusters K exceeds 
the number of clusters present in the data. Following the concept of sparse finite mixtures 
(Malsiner-Walli et ah, 2016), we choose a symmetric Dirichlet distribution as prior for 
the weight distribution, i.e. rj\eo ~ Dirx(e o), and base our choice of eo on the results 
of Rousseau and Mengersen (2011) concerning the asymptotic behavior of the posterior 
distribution of an overfitting mixture model. They show that this behavior is determined 
by the hyperparameter eo of the Dirichlet prior on the weights. In particular, they prove 
that, if eo < d/2, where d is the dimension of the cluster-specific parameters Ok, then the 
posterior expectation of the weights associated with superfluous clusters asymptotically 
converges to zero. 

Hence, we specify a sparse prior on the cluster weights 77 by choosing eo <C d/2 so that 
superfluous clusters are emptied during MCMC sampling and the number of non-empty 
clusters on the cluster level is an estimator for the unknown number of data clusters. In 
this way, the specification of a sparse cluster weight prior in an overfitting mixture of 
mixtures model provides an “automatic tool” to select the number of clusters, avoiding 
the expensive computation of marginal likelihoods as, e.g., in Friihwirth-Schnatter (2004). 
Empirical results in Malsiner-Walli et al. (2016) indicate that eo needs to be chosen very 
small, e.g. e 0 = 0 . 001 , to actually empty all superfluous clusters in the finite sample case. 

On the lower level (2), in each cluster k, a semi-parametric approximation of the 
cluster distributions is achieved by mixing L multivariate Gaussian subcomponent den¬ 
sities /V(y;!/•%, ^ki), l = 1 according to the subcomponent weight vector w fc = 
{wk 1 , • • •, t Uki/)- The cluster-specihc parameter vector 

Ok = ( w fci O'kli ■ ■ ■ ) O'kLi • ■ • 1 ^hli) (3) 

consists of Wfc as well as the means 7 % and covariance matrices of all Gaussian sub¬ 
component densities. L is typically unknown, but as we are not interested in estimating the 
“true” number of subcomponents L forming the cluster, we only ensure that L is chosen 
sufficiently large to obtain an accurate approximation of the cluster distributions. While 
the choice of L is not crucial to ensure a good model fit as long as L is sufficiently large, a 
too generous choice of L should be avoided for computational reasons as the computational 
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complexity of the estimation increases with the number of subcomponents L. 

By choosing the prior ~ Dir^do ) with do = d/2+2, the approximation of the cluster 
density is obtained by filling all L subcomponents, thus avoiding empty subcomponents. 
This choice is motivated again by the results of Rousseau and Mengersen (2011) who show 
that, if do > d/2, the posterior density asymptotically handles an overfitting mixture by 
splitting “true” components into two or more identical components. 

2.2 Identification through hierarchical priors 

When fitting the finite mixture model (1) with semi-parametric cluster densities given 
by (2), we face a special identifiability problem, since the likelihood is entirely agnostic 
about which subcomponents form a cluster. Indeed, the likelihood is completely ignorant 
concerning the issue which of the K-L components belong together, since (1) can be written 
as an expanded Gaussian mixture with K-L components with weights w k i = T)k w ki, 

K L 

p(y*l©>*7) = EE w k if+{yi\^ki^ki)- (4) 

k =1 1=1 

These K-L components can be permuted in ( K-L) \ different ways and the resulting ordering 
can be used to group them into K different cluster densities, without changing the mixture 
likelihood (4). Hence, the identification of (1), up to label switching on the upper level, 
hinges entirely on the prior distribution. 

Subsequently, we suggest a hierarchical prior that addresses these issues explicitly. Con¬ 
ditional on a set of fixed hyperparameters (p 0 = (e 0 , d 0 , c 0 , go, G 0 , B 0 , m 0 , M 0 , u), the weight 
distribution rj\eo ~ Dirx{e o) and the K cluster-specific parameter vectors Qk\4>o ~ P(0 k I0o) 
are independent a priori, i.e.: 

K 

p(rj,e 1 ,... ,e K \(j)o) = p(r]\eo)Y[p(O k \(j)o). (5) 

k =1 

This prior formulation ensures that the K non-Gaussian cluster distributions of the upper 
level mixture (1) are invariant to permutations. Within each cluster k, the prior distribution 
p(O k \</>o) admits the following block independence structure: 

p(0fc|0o) = p(w k \d 0 )p{ii kl ,.. •, MfciJ B o, m o, M 0 , z'MEfci,..., T, kL \ c 0 , g 0 , G 0 ), (6) 
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where Wfe|d 0 ~ DirL(d 0 ). Conditional on 0 O , the subcomponent means /x fel ,... , n kL are de¬ 
pendent a priori as are the subcomponent covariance matrices ..., However, they 
are assumed to be exchangeable to guarantee that within each cluster k, the L Gaussian 
subcomponents in (2) can be permuted without changing the prior. 

To create this dependence, a hierarchical “random effects” prior is formulated, where, on 
the upper level, conditional on the fixed upper level hyperparameters (g 0 , G 0 , m 0 , M 0 , z/), 
cluster specific random hyperparameters (Cofc, b 0 fc), and A k = diag(Afci,..., Afc r ), are drawn 
independently for each k — 1,..., K from a set of three independent base distributions: 

Cofc|<?o, G 0 ~ W r (< 7 o, G 0 ), b 0 fc|m 0 , M 0 ~ A/" r (m 0 , M 0 ), (A fel ,..., A fcr )| v~g(v,v), (7) 

where Af r () and W r () denote the r-multivariate normal and Wishart distribution, respec¬ 
tively, and Q() the gamma distribution, parametrized such that E(\ki\v) = 1. 

On the lower level, conditional on the cluster specific random hyperparameters (Cofc, b 0 fc, 
Afe) and the fixed lower level hyperparameters (B 0 , Co), the L subcomponent means n kl and 
covariance matrices are drawn independently for all l — 1,..., L: 

Mfc/1B 0 , b 0 fc, A k ~ A/" r (b ofc , y / A^B 0 y / A/?), T, kl 1 \c 0: C ok ~ hV r (c 0 , Co*,). (8) 

2.3 Tuning the hyperparameters 

To identify the mixture of mixtures model given in (1) and (2) through the prior defined 
in Section 2.2, the fixed hyperparameters 0o have to be chosen carefully. In addition, we 
select them in a way to take the data scaling into account, avoiding the need to standardize 
the data prior to data analysis. 

First, it is essential to clarify what kind of shapes and forms are aimed at as cluster 
distributions. We give the following (vague) characterization of a data cluster: A data 
cluster is a very “dense” region of data points, with possibly no “gaps” within the cluster 
distribution, whereas different clusters should be located well-separated from each other, 
i.e. here large “gaps” between the cluster distributions are desired. We confine ourselves to 
the investigation of clusters with approximately convex cluster shapes, where the cluster 
center can be seen as a suitable representative for the entire cluster. Regarding volume, 


8 


orientation or asymmetry of the data clusters we are looking for, no constraints on the 
cluster shapes and forms are imposed. 

Based on this cluster concept, our aim is to model a dense and connected cluster distri¬ 
bution by a mixture of normal subcomponents. Various strategies regarding the modeling 
of the subcomponent means and covariance matrices could be employed. We decided to 
allow for flexible shapes for the single subcomponents, ensuring that they strongly overlap 
at the same time. An alternative approach would be to use constrained simple shaped 
subcomponents, e.g., subcomponents with isotropic covariance matrices. However, in this 
case a large number of subcomponents might be needed to cover the whole cluster region 
and shrinkage of the subcomponent means toward the common cluster center may not be 
possible. Since then some of the subcomponents have to be located far away from the clus¬ 
ter center in order to fit also boundary points, considerable distances have to be allowed 
between subcomponent means. This induces the risk of gaps within the cluster distribution 
and a connected cluster distribution may not result. Therefore, in our approach the cluster 
distributions are estimated as mixtures of only a few but unconstrained, highly dispersed 
and heavily overlapping subcomponents where the means are strongly pulled toward the 
cluster center. In this way, a connected cluster distribution is ensured. 

In a Bayesian framework, we need to translate these modeling purposes into appropriate 
choices of hyperparameters. On the upper level, the covariance matrix M 0 controls the 
amount of prior shrinkage of the cluster centers bofc toward the overall data center mo, 
which we specify as the midpoint of the data. To obtain a prior, where the cluster centers 
b 0 fc are allowed to be widely spread apart and almost no shrinkage toward m 0 takes place, 
we choose M 0 S y , where S y is the sample covariance matrix of all data, e.g. M 0 = 10S y . 

Our strategy for appropriately specifying the hyperparameters Go and B 0 is based on 
the variance decomposition of the mixture of mixtures model, which splits Cov (Y) into the 
different sources of variation. For a finite mixture model with K clusters, as given in (1), the 
total heterogeneity Cov(Y) can be decomposed in the following way (Friihwirth-Schnatter, 
2006, p. 170): 

K K 

Cov(Y) = ^ Tj k n k p,' k - /i/V = (1 - <f> B )Cov{ Y) + (p B Cov( Y), (9) 

k= 1 k =1 

where the cluster means n k and the cluster covariance matrices X*, are the first and second 
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moments of the cluster distribution Pk{yi\6k) and /i = Ylk r)kH k is the mixture mean. In this 
decomposition (f> B is the proportion of the total heterogeneity explained by the variability 
of the cluster means and (1 — (p B ) is the proportion explained by the average variability 
within the clusters. The larger (f>B, the more the clusters are separated, as illustrated in 
Figure 1 for a three-component standard Gaussian mixture with varying values of <p B ■ 

For a mixture of mixtures model, the heterogeneity (1 — (p B )Cov(Y) explained within 
a cluster can be split further into two sources of variability, namely the proportion (f>w 
explained by the variability of the subcomponent means fi kl around the cluster center fi k , 
and the proportion (1 ~4>w) explained by the average variability within the subcomponents: 

K K 

Cov(Y) = E rjk'Ek + Y - aX 

k= 1 k= 1 

K L K / L \ K 

= 5Z ^ W ^kl + Y Tlk ( Y Wk ^kl^'kl ~ MiX I + Y WkK ~ W' (10) 

k =1 1=1 k =1 \ /=1 / k =1 

= (1 - <j, w ){ 1 - Y) + 4> w { 1 - <p B )Cov{ Y) + <p B Cov{ Y). 

Based on this variance decomposition we select the proportions <p B and (f>w and incorporate 
them into the specification of the hyperparameters of our hierarchical prior. 

(f) B defines the proportion of variability explained by the different cluster means. We 
suggest to specify <p B not too large, e.g., to use 4> B = 0.5. This specification may seem to 
be counterintuitive as in order to model well-separated clusters it would seem appropriate 
to select (ft B large. However, if (f> B is large, the major part of the total heterogeneity of the 
data is already explained by the variation (and separation) of the cluster means, and, as a 
consequence, only a small amount of heterogeneity is left for the within-cluster variability. 
This within-cluster variability in turn will get even more diminished by the variability 
explained by the subcomponent means leading to a small amount of variability left for the 
subcomponents. Thus for large values of <p B , estimation of tight subcomponent densities 
would result, undermining our modeling aims. 

(pw defines the proportion of within-cluster variability explained by the subcomponent 
means, (bw also controls how strongly the subcomponent means are pulled together and 
influences the overlap of the subcomponent densities. To achieve strong shrinkage of the 
subcomponent means toward the cluster center, we select small values of <pw, e.g. 4>w = 0.1. 
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Figure 1: Variance decomposition of a mixture distribution. Scatter plots of samples from a standard 
normal mixture distribution with three components and equal weights, with a varying amount of hetero¬ 
geneity 4>b explained by the variation of the component means, </>b = 0.1, 4>b = 0.5 and (f>B = 0.9 (from 
left to right). 

Larger values of (pw may introduce gaps within a cluster, which we want to avoid. 

Given <pB and (pw , we specify the scale matrix Go of the prior on Cofc such that the 
a priori expectation of the first term in the variance decomposition (10), given by 

( k l \ K L 

j = '^E(ri k )^^E(w k i)E(E{Y lk i\C ok )) = ^ 0 /(c 0 - (r + 1)/2)G 0 " 1 , 

k =1 1=1 J k = 1 1=1 

matches the desired amount of heterogeneity explained by a subcomponent: 

W(co - {r + 1)/2)Gq 1 = (1 - <p w )( 1 - <p B )Cov( Y). (11) 

We replace Cov(Y) in (11) with the main diagonal of the sample covariance S y to take 
only the scaling of the data into account (see e.g. Fruhwirth-Schnatter, 2006). This gives 
the following specification for Go: 

Gq 1 = (1 - 0w)(l - 0 b)(c o - (r + l)/2)/g 0 ■ diag(S„). (12) 

Specification of the prior of the subcomponent covariance matrices S k i, ■ ■ ■, S k L is com¬ 
pleted by defining the scalar prior hyperparameters Co and r/ 0 ■ Fruhwirth-Schnatter (2006, 
Section 6.3.2, p. 192) suggests to set c 0 > 2 + (r — l)/2. In this way the eigenvalues of 
'Ski'^km are b° un( ied away from 0 avoiding singular matrices. We set Co = 2.5 + (r — l)/2 
to allow for a large variability of The Wishart density is regular if go > (r — l)/2 and 
in the following we set go = 0.5 + (r — l)/2. 
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Regarding the prior specification of the subcomponent means fi kl , ..., n kL , we select the 
scale matrix Bq in order to concentrate a lot of mass near the cluster center bofc, pulling 
fi kl toward bofc. Matching the a priori expectation of the second term in the variance 
decomposition (10), given by 

E I rik ( ^2 w ^ki^'ki ~ AR/4 = ^2 E (Vk)^2 E(w k i)E(^ kl n' kl - n k /j,' k ) = B 0 , 

\k =1 \l =1 // k =1 1=1 

to the desired proportion of explained heterogeneity and, using once more only the main 
diagonal of S y we obtain B 0 = — 0s)diag(S y ), which incorporates our idea that only a 

small proportion <pw of the within-cluster variability should be explained by the variability 
of the subcomponent means. 

After having chosen (ft B and basically the cluster structure and shape is a priori 
determined. However, in order to allow for more flexibility in capturing the unknown cluster 
shapes in the sense that within each cluster the amount of shrinkage of the subcomponent 
means fx kl toward the cluster center bofc need not to be the same for all dimensions, for each 
cluster k and each dimension j additionally a random adaptation factor A k j is introduced 
in (8) which adjusts B 0 . The gamma prior for A k j in (7) implies that the prior expectation 
of the covariance matrix of n kl equals B 0 . However, A k j acts as a local adjustment factor 
for cluster k which allows to shrink (or inflate) the variance of subcomponent means fi k ij in 
dimension j in order to adapt to a more (or less) dense cluster distribution as specified by 
Bo- In order to allow only for small adjustments of the specified Bo, we choose zy = 10, in 
this way almost 90% of the a priori values of A k j are between 0.5 and 1.5. This hierarchical 
prior specification for n kl corresponds to the normal gamma prior (Griffin and Brown, 2010) 
which has been applied by Friihwirth-Schnatter (2011) and Malsiner-Walli et al. (2016) in 
the context of finite mixture models for variable selection. 

2.4 Relation to BNP mixtures 

Our approach bears resemblance to various approaches in BNP modeling. First of all, the 
concept of sparse finite mixtures as used in Malsiner-Walli et al. (2016) is related to Dirichlet 
process (DP) mixtures (Muller and Mitra, 2013) where the discrete mixing distribution in 
the finite mixture (1) is substituted by a random distribution G ~ DP(a, H ), drawn from 
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a DP prior with precision parameter a and base measure H. As a draw G from a DP is 
almost surely discrete, the corresponding model has a representation as an infinite mixture: 

OO 

p(y) = ^2vkPk{y\0k), (13) 

fc=l 

with i.i.d. atoms Gk ~ H drawn from the base measure H and weights r/k = Vk n >=i (1 — v j) 
obeying the stick breaking representation with ~ B (l,a) (Sethuraman, 1994). 

If the hyperparameter in the weight distribution r] of a sparse finite mixture is chosen 
as eo = a/K, i.e. rj ~ Dirx{oi/K), and the component parameters Gk ~ H are i.i.d. 
draws from H, then as K increases, the sparse finite mixture in Equation (1) converges to 
a DP mixture with mixing distribution G ~ DP(a, H), see Green and Richardson (2001). 
For example, the sparse finite Gaussian mixture introduced in Malsiner-Walli et al. (2016) 
converges to a Dirichlct process Gaussian mixture as K increases, with (fi k , 5R ; ) being i.i.d. 
draws from the appropriate base measure H. 

The more general sparse finite mixture of mixtures model introduced in this paper also 
converges to a Dirichlet process mixture where the atoms are finite mixtures indexed by 
the parameter Gk defined in (3). The parameters Gk are i.i.d. draws from the base measure 
(6), with strong dependence among the means /j , kl ,..., fi kL and covariances S fcl ,..., H k L 
within each cluster k. This dependence is achieved through the two-layer hierarchical 
prior described in (7) and (8) and is essential to create well-connected clusters from the 
subcomponents, as outlined in Section 2.3. 

Also in the BNP framework models have been introduced that create dependence, either 
in the atoms and/or in the weights attached to the atoms. For instance, the nested DP 
process of Rodriguez et al. (2008) allows to cluster distributions across N units. Within 
each unit i, i = 1, repeated (univariate) measurements y it ,t = 1 ,...,W arise as 

independent realizations of a DP Gaussian mixture with random mixing distribution G*. 
The GiS are i.i.d. draws from a DP, in which the base measure is itself a Dirichlet process 
DP(/3, H), i.e. Gi ~ DP(a , DP(/3, H)). Hence, two distributions Gi and Gj either share the 
same weights and atoms sampled from H , or the weights and atoms are entirely different. 
If only a single observation r/j is available in each unit, i.e. N, = 1, then the nested DP 
is related to our model. In particular, it has a two-layer representation as in (1) and (2), 
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however with both K and L being infinite. The nested DP can, in principal, be extended 
to multivariate observations y*. In this case, p(yi) takes the same form as in (13), with the 
same stick breaking representation for the cluster weights 771 , 772 , ■ • •• On the lower level, 
each cluster distribution p k (yi\6 k ) is a DP Gaussian mixture: 

OO 

Pk(yi\Ok) = ^WkifAfiyilVkh s fcz), (14) 

1=1 

where the component weights w k i are derived from the stick breaking representation w k i = 
u k i n^=i(! ~ u kj), 1 = 1,2,... where u k] 13(1,13). For the nested DP, dependence is intro¬ 
duced only on the level of the weights and sticks, as the component parameters pi kl , ~ H 

are i.i.d. draws from the base measure H. This lack of prior dependence among the atoms 
(H k i, Sfc/) is likely to be an obstacle in a clustering context. 

The BNP approach most closely related to our model is the infinite mixture of infinite 
Gaussian mixtures (I 2 GMM) model of Yerebakan et al. (2014) which also deals with clus¬ 
tering multivariate observations from non-Gaussian component densities. 1 The fiGMM 
model has a two-layer hierarchical representation like the nested DP. On the top level, i.i.d. 
cluster specific locations b 0 A, and covariances S*. are drawn from a random distribution 
G ~ I)P ( a . H) arising from a DP prior with base measure H being equal to the conju¬ 
gate normal-inverse-Wishart distribution. A cluster specific DP is introduced on the lower 
level as for the nested DP; however, the I 2 GMM model is more flexible, as prior depen¬ 
dence is also introduced among the atoms belonging to the same cluster. More precisely, 
y j ~ A f r (lii, S fc ), with /x.- ~ G k , where G k ~ DP((3, H k ) is a draw from a DP with cluster 
specific base measure H k = Af r (h ok , 1 ). 

It is easy to show that the I 2 GMM model has an infinite two-layer representation as in 
(13) and (14), with exactly the same stick breaking representation. 2 However, the I 2 GMM 
model has a constrained form on the lower level, with homoscedastic covariances = Sfc, 
whereas the locations fi kl scatter around the cluster centers b 0 fc as in our model: 

(b ofc , £*) ~ H, /xjbofc, S, ~ H k . (15) 

would like to thank a reviewer for pointing us to this paper. 

2 Note that the notation in Yerebakan et al. (2014) is slightly different, with 7 and a corresponding to 

a and (3 introduced above. 
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In our sparse mixture of mixtures model, we found it useful to base the density estimator on 
heteroscedastic covariances 'Em, to better accommodate the non-Gaussianity of the cluster 
densities with a fairly small number L of subcomponents. It should be noted that our 
semi-parametric density estimator is allowed to display non-convex shapes, as illustrated 
in Figure C.2 in the Appendix. Nevertheless, we could have considered a mixture in (2) 
where £&/ = £&, with the same base measure for the atoms (/tx fcl ,..., fJ. kL , £&) as in (15). 
In this case, the relationship between our sparse finite mixture and the I 2 GMM model 
would become even more apparent: by choosing e 0 = a/K and d 0 = fd/L and letting K 
and L go to infinity, our model would converge to the I 2 GMM model. 

3 Clustering and posterior inference 

3.1 Clustering and selecting the number of clusters 

For posterior inference, two sequences of allocation variables are introduced, namely the 
cluster assignment indicators S = (Si,... ,Sn) and the within-cluster allocation variables 
I = (/ 1 ,..., /at). More specifically, S) G {1,..., K } assigns each observation y,; to cluster 
Si on the upper level of the mixture of mixtures model. On the lower level, G {1,..., L} 
assigns observation y ?: to subcomponent Hence, the pair (Si, Ij) carries all the informa¬ 
tion needed to assign each observation to a unique component in the expanded mixture 

(4). 

Note that for all observations y, and y 3 belonging to the same cluster, the upper level 
indicators S) = Sj will be the same, while the lower level indicators I L ^ Ij might be 
different, meaning that they belong to different subcomponents within the same cluster. 
It should be noted that the Dirichlet prior w k ~ Diri(d 0 ), with d 0 > d/2, on the weight 
distribution ensures overlapping densities within each cluster, in particular if L is overfit¬ 
ting. Hence the indicators A will typically cover all possible values {1,... ,L} within each 
cluster. 

For clustering, only the upper level indicators S are explored, integrating implicitly 
over the uncertainty of assignment to the subcomponents on the lower level. A cluster 
C k = {i\Si = k} is thus a subset of the data indices {!,..., N}, containing all observations 
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with identical upper level indicators. Hence, the indicators S define a random partition 
V = {Ci,..., Cr 0 } of the N data points in the sense of Lau and Green (2007), as y \ and 
y j belong to the same cluster, if and only if S t = Sj. The partition V contains AT 0 = \P\ 
clusters, where \V\ is the cardinality of V. Due to the Dirichlet prior 77 ~ Dir^^e 0 ), with 
eo close to 0 to obtain a sparse finite mixture, K 0 is a random number being a priori much 
smaller than K. 

For a sparse finite mixture model with K clusters, the prior distribution over all ran¬ 
dom partitions V of N observations is derived from the joint (marginal) prior p( S) = 
/ nf= 1 p(S i \r})dr] which is given, e.g., in Friihwirth-Schnatter (2006, p. 66): 


P( s) 


r(AT 0 ) 

r(N + Ke 0 )r(e 0 ) K ° 


J] r(Nk + e 0 ), 

k:N k > 0 


(16) 


where iV* = #{Si = k}. For a given partition V with K 0 data clusters, there are K\/{K — 
A'o)! assignment vectors S that belong to the equivalence class defined by V. The prior 
distribution over all random partitions V is then obtained by summing over all assignment 
vectors S that belong to the equivalence class defined by V: 

?(P|A ' o) = (K - K 0 )\ UN + L„)r(e„)& fc Il >o r (^ + e »>- ( 17 > 

which takes the form of a product partition model and therefore is invariant to permuting 
the cluster labels. Hence, it is possible to derive the prior predictive distribution p(*S' i |S_*), 
where S _j denote all indicators, excluding Si. Let Kq 1 be the number of non-empty clusters 
implied by S_j and let N^ 1 be the corresponding cluster sizes. From (16), we obtain the 
following probability that Si is assigned to an existing cluster k\ 

Pr{S s = *|S-„ 1 V > 0} = (18) 

The prior probability that Si creates a new cluster with Si G / = {k\N^ 1 = 0} is equal to 

PrG 7|S_i} = (K - K^)Pr{Si = A;*|S_<, k* G 1} = C ° (A ~ A °!l (19) 

iv — i + eoA 


It is illuminating to investigate the prior probability to create new clusters in detail. First 
of all, for e 0 independent of K , this probability not only depends on e 0 , but also increases 
with K. Hence a sparse finite mixture model based on the prior 77 ~ T>k (eo) can be 


16 






regarded as a two-parameter model, where both e 0 and K influence the a priori expected 
number of data clusters Kq which is determined for a DP mixture solely by a. A BNP 
two-parameter mixture is obtained from the Pitman-Yor process (PYP) prior PY((3, a ) 
with f3 G [0,1), a > —(3 (Pitman and Yor, 1997), with stickbreaking representation Vk ~ 
B (1 — /3, a + k/3). The DP prior results as that special case where (3 = 0. 

Second, the prior probability (19) to create new clusters in a sparse finite mixture model 
decreases, as the number Kq 1 of non-empty clusters increases. This is in sharp contrast to 
DP mixtures where this probability is constant and PYP mixtures where this probability 
increases, see e.g., Fall and Barat (2014). 

Finally, what distinguishes a sparse finite mixture model, both from a DP as well as 
a PYP mixture, is the a priori expected number of data clusters K 0 , as the number N 
of observations increases. For K and eo independent of N, the probability to create new 
clusters decreases, as N increases, and converges to 0, as N goes to infinity. Therefore, K 0 
is asymptotically independent of N for sparse finite mixtures, whereas for the DP process 
K 0 ~ a log(iV) (Korwar and Hollander, 1973) and K 0 ~ N 13 obeys a power law for PYP 
mixtures (Fall and Barat, 2014). This leads to quite different clustering behavior for these 
three types of mixtures. 

A well-known limitation of DP priors is that a priori the cluster sizes are expected to 
be geometrically ordered, with one big cluster, geometrically smaller clusters, and many 
singleton clusters (Muller and Mitra, 2013). PYP mixtures are known to be more useful 
than the DP mixture for data with many significant, but small clusters. A common criticism 
concerning finite mixtures is that the number of clusters needs to be known a priori. Since 
this is not the case for sparse finite mixtures, they are useful in the context of clustering, 
in particular in cases where the data arise from a moderate number of clusters, that does 
not increase as the number of data points N increases. 

3.2 MCMC estimation and posterior inference 

Bayesian estimation of the sparse hierarchical mixture of mixtures model is performed using 
MCMC methods based on data augmentation and Gibbs sampling. We only need standard 
Gibbs sampling steps, see the detailed MCMC sampling scheme in Appendix A. 
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In order to perform inference based on the MCMC draws, i.e. to cluster the data, to 
estimate the number of clusters, to solve the label switching problem on the higher level and 
to estimate cluster-specific parameters, several existing procedures can be easily adapted 
and applied to post-process the posterior draws of a mixture of mixtures model, e.g., those 
which are, for instance, implemented in the R packages PReMiuM (Liverani et al., 2015) 
and labehswitching (Papastamoulis, 2015). 

For instance, the approach in PReMiuM is based on the posterior probabilities of co¬ 
clustering, expressed through the similarity matrix PrlS 1 * = Sjly} which can be estimated 
from the M posterior draws S= 1,..., M , see Appendix B for details. The methods 
implemented in labehswitching aim at resolving the label switching problem when fitting 
a finite mixture model using Bayesian estimation. Note that in the case of the mixture 
of mixtures model label switching occurs on two levels. On the cluster level, the label 
switching problem is caused by invariance of the mixture likelihood given in Equation (1) 
with respect to reordering of the clusters. On this level, label switching has to be resolved, 
since the single cluster distributions need to be identified. On the subcomponent level, label 
switching happens due to the invariance of Equation (2) with respect to reordering of the 
subcomponents. As we are only interested in estimating the entire cluster distributions, 
it is not necessary to identify the single subcomponents. Therefore, the label switching 
problem can be ignored on this level. 

In this paper, the post-processing approach employed first performs a model selection 
step. The posterior draws of the indicators S ( m \m — 1,..., M are used to infer the number 
of non-empty clusters on the upper level of the mixture of mixtures model and the 
number of data clusters is then estimated as the mode. Conditional on the selected model, 
an identified model is obtained based on the point process representation of the estimated 
mixture. This method was introduced in Friihwirth-Schnatter (2006, p. 96) and successfully 
applied to model-based clustering in various applied research, see e.g. Friihwirth-Schnatter 
(2011) for some review. This procedure has been adapted to sparse finite mixtures in 
Friihwirth-Schnatter (2011) and Malsiner-Walli et al. (2016) and is easily extended to deal 
with sparse mixture of mixtures models, see Appendix B for more details. We will use 
this post-processing approach in our simulation studies and the applications in Section 4 


18 


and Appendices C, D and F to determine a partition of the data based on the maximum 
a posteriori (MAP) estimates of the relabeled cluster assignments. 

4 Simulation studies and applications 

The performance of the proposed strategy for selecting the unknown number of clusters 
and identifying the cluster distributions is illustrated in two simulation studies. In the 
first simulation study we investigate whether we are able to capture dense non-Gaussian 
data clusters and estimate the true number of data clusters. Furthermore, the influence 
of the specified maximum number of clusters K and subcomponents L on the clustering 
results is studied. In the second simulation study the sensitivity of the a priori defined 
proportions <f>B and (ftw on the clustering result is investigated. For a detailed description 
of the simulation design and results see Appendix C. Overall, the results indicated that our 
approach performed well and yielded promising results. 

To further evaluate our approach, we fit the sparse hierarchical mixture of mixtures 
model on benchmark data sets and real data. First, we consider five data sets which 
were previously used to benchmark algorithms in cluster analysis. For these data sets we 
additionally apply the “merging strategy” proposed by Baudry et al. (2010) in order to 
compare the results to those of our approach. For these benchmark data sets class labels 
are available and we assess the performance by comparing how well our approach is able 
to predict the class labels using the cluster assignments, measured by the misclassihcation 
rate as well as the adjusted Rand index. 

To assess how the algorithm scales to larger data sets we investigate the application to 
two flow cytometry data sets. The three-dimensional DLBCL data set (Lee and McLachlan, 
2013b) consists of around 8000 observations and comes with manual class labels which can 
be used as benchmark. The GvHD data set (Brinkman et ah, 2007) consists of 12441 
observations, but no class labels are available. We compare the clusters detected for this 
data set qualitatively to solutions previously reported in the literature. 

The detailed description of all investigated data sets as well as of the derivation of the 
performance measures are given in Appendix D. For the benchmark data sets, the number 
of estimated clusters A'o, the adjusted Rand index ( adj ), and misclassihcation rate (er) 
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are reported in Table 1 for all estimated models. In the first columns of Table 1, the name 
of the data set, the number of observations N, the number of variables r and the number 
of true classes K true (if known) are reported. To compare our approach to the merging 
approach proposed by Baudry et al. (2010), we use the function Mclust of the R package 
mclust (Fraley et ah, 2012) to first fit a standard normal mixture distribution with the 
maximum number of components K = 10. The number of estimated normal components 
based on the BIC is reported in the column Mclust. Then the selected components are 
combined hierarchically to clusters by calling function clustCombi from the same package 
(column clustCombi). The number of clusters is chosen by visual detection of the change 
point in the plot of the rescaled differences between successive entropy values, as suggested 
by Baudry et al. (2010). Furthermore, to compare our results to those obtained if a cluster 
distribution is modeled by a single normal distribution only, a sparse finite mixture model 
with K = 10 (Malsiner-Walli et ah, 2016) is fitted to the data sets (column SparseMix ). 
The results of fitting a sparse hierarchical mixture of mixtures model with K = 10 are 
given in column Spars eMixMix, where L = 5 is compared to our default choice of L = 4 
to investigate robustness with respect to the choice of L. For each estimation, MCMC 
sampling is run for 4000 iterations after a burn-in of 4000 iterations. 

As can be seen in Table 1, for all data sets the sparse hierarchical mixture of mixtures 
model is able to capture the data clusters quite well both in terms of the estimated number 
of clusters and the clustering quality measured by the misclassihcation rate as well as the 
adjusted Rand index. In general, our approach is not only outperforming the standard 
model-based clustering model using mixtures of Gaussians regarding both measures, but 
also the approach proposed by Baudry et al. (2010). In addition, it can be noted that for 
all data sets the estimation results remain quite stable, if the number of subcomponents 
L is increased to 5, see the last column in Table 1. The results for the Yeast data set 
are of particular interest as they indicate that clustCombi completely fails. Although the 
misclassihcation rate of 25% implies that only a quarter of the observations is assigned to 
“wrong” clusters, inspection of the clustering obtained reveals that almost all observations 
are lumped together in a single, very large cluster, whereas the few remaining observations 
are split into five very small clusters. This bad clustering quality is better reflected by the 
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Mclust 

SparseMix 

SparseMixMix 





K 

= 10 

I\ = 10 


K = 

10 

Data set 

N 

r 

j^true 

Mclust 

clustCombi 

L = 1 

L = 4 


L = 5 

Yeast 

626 

3 

2 

8 (.50, .20) 

6 (-.02, 0.25) 

6 (48, .23) 

2 (.68, 

08) 

2 (.71, .07) 

Flea beetles 

74 

6 

3 

5 (.77, .18) 

4 (.97, .03) 

3 (1.00, .00) 

3 (1.00, 

.00) 

3 (1, .00) 

AIS 

202 

3 

2 

3 (.73, .13) 

2 (.66, .09) 

3 (.76, .11) 

2 (.81, 

05) 

2 (.76, .06) 

Wisconsin 

569 

3 

2 

4 (.55, .30) 

4 (.55, .30) 

4 (.62, .21) 

2 (.82, 

05) 

2 (.82, .05) 

Flower 

400 

2 

4 

6 (.52, .35) 

4 (.99, .01) 

5 (.67, .20) 

4 (.97, 

01) 

4 (.97, .02) 


Table 1: Results for the estimated number of data clusters K$ for various benchmark data sets, using 
the functions Mclust to fit a standard mixture model with K = 10 and clustCombi to estimate a mixture 
with combined components (column Mclust ), using a sparse finite mixture model with K = 10 (column 
SparseMix), and estimating a sparse hierarchical mixture of mixtures model with K = 10, (j)B = 0.5 and 
<t>w = 0.1, and L = 4, 5 (column SparseMixMix). Priors and hyperparameter specifications are selected as 
described in Section 2. In parentheses, the adjusted Rand index (“1” corresponds to perfect classification) 
and the proportion of misclassified observations (“ 0 ” corresponds to perfect classification) are reported. 

adjusted Rand index which takes a negative value ( adj = —0.02), i.e. is “worse than would 
be expected by guessing” (Franczak et ah, 2012). For the flower data set, more results are 
given in Appendix D where the obtained clustering and cluster distributions are illustrated. 

In order to investigate the performance of our approach on larger data sets with a 
slightly different cluster structure, we fit the sparse hierarchical mixture of mixtures model 
to two flow cytometry data sets. These applications also allow us to indicate how the 
prior settings need to be adapted if a different cluster structure is assumed to be present 
in the data. As generally known, flow cytometry data exhibit non-Gaussian characteristics 
such as skewness, multimodality and a large number of outliers, as can be seen in the 
scatter plot of two variables of the GvHD data set in Figure 3. Thus, we specified a sparse 
hierarchical mixture of mixtures model with K = 30 clusters and increased the number of 
subcomponents forming a cluster to L — 15 in order to handle more complex shapes of 
the cluster distributions given the large amount of data. Since the flow cytometry data 
clusters have a lot of outliers similar to the clusters generated by shifted asymmetric Laplace 
(SAL) distributions (see Appendix F), we substitute the hyperprior C 0 fc ~ VV r (r/ 0 , G 0 ) by 
the fixed value Cofc = ^oGq 1 and set A kj = 1 , j — 1 ,... , r to prevent that within a cluster 
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Figure 2: Flow cytometry data set DLBCL. Scatterplot of the clustering results. 

the subcomponent covariance matrices are overly shrunken and become too similar. In this 
way, subcomponent covariance matrices are allowed to vary considerably within a cluster 
and capture both a dense cluster region around the cluster center and scattered regions at 
the boundary of the cluster. 

We fit this sparse hierarchical mixture of mixtures model to the DLBCL data after 
removing 251 dead cells. For most MCMC runs after a few hundred iterations all but 
four clusters become empty during MCMC sampling. The estimated four cluster solution 
coincides almost exactly with the cluster solution obtained with manual gating; the adjusted 
Rand index is 0.95 and the error rate equals 0.03. This error rate outperforms the error rate 
of 0.056 reported by Lee and McLachlan (2013b). In Figure 2 the estimated four cluster 
solution is visualized. 

When fitting a sparse hierarchical mixture of mixtures model to the GvHD data, the 
classifications resulting from different runs of the MCMC algorithm seemed to be rather 
stable. The obtained solutions differ mainly in the size of the two large clusters with low 
expressions. These, however, are supposed to not contain any information regarding the 
development of the disease. On the right hand side of Figure 3, the results of one specific 
run are shown in a heatmap. In this run, we found eight clusters which are similar to those 
reported by Friihwirth-Schnatter and Pyne (2010) when fitting a skew-f mixture model to 
these data. In the heatmap each row represents the location of a six-dimensional cluster, 
and each column represents a particular marker (variable). The red, white and blue colors 
denote high, medium and low expressions. 

As in Friihwirth-Schnatter and Pyne (2010), we identified two larger clusters (43% and 
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Figure 3: Flow cytometry data set GvHD. Scatter plot of two variables (“FSC”, “CD8”) (left-hand side), 
and heatmap of the clustering results by fitting a sparse hierarchical mixture of mixtures model (right- 
hand side). In the heatmap, each row represents the location of a six-dimensional cluster, and each column 
represents a particular marker. The red, white and blue colors denote high, medium and low expression, 
respectively. 


20.4%, first two rows in the heatmap) with rather low expressions in the last four variables. 
We also identified a smaller cluster (3.8%, forth row from the bottom) representing live cells 
(high values in the first two variables) with a unique signature in the other four variables 
(high values in all four variables). Also two other small clusters can be identified (second 
and third row from the bottom) which have a signature very similar to the clusters found 
by Friihwirth-Schnatter and Pyne (2010), and thus our results confirm their findings. 


5 Discussion 

We propose suitable priors for fitting an identified mixture of normal mixtures model within 
the Bayesian framework of model-based clustering. This approach allows for (1) automatic 
determination of the number of clusters and (2) semi-parametric approximation of non- 
Gaussian cluster distributions by mixtures of normals. We only require the assumption that 
the cluster distributions are dense and connected. Our approach consists in the specification 
of structured informative priors on all model parameters. This imposes a rigid hierarchical 
structure on the normal subcomponents and allows for simultaneous estimation of the 
number of clusters and their approximating distributions. This is in contrast to the two- 
step merging approaches, where in the first step the data distribution is approximated by 
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a suitable normal mixture model. However, because this approximation is made without 
taking the data clusters into account which are reconstructed only in the second step of 
the procedure, the general cluster structure might be missed by these approaches. 

As we noted in onr simulation studies, the way in which the cluster mixture distributions 
are modeled by the subcomponent densities is crucial for the clustering result. Enforcing 
overlapping subcomponent densities is essential in order to avoid that a single subcom¬ 
ponent becomes too narrow thus leading to a small a posteriori cluster probability for 
observations from this subcomponent. Also, enforcing that observations are assigned to all 
subcomponents during MCMC sampling is important as the estimation of empty subcom¬ 
ponents would bias the resulting cluster distribution because of the “prior” subcomponents. 
For modeling large, overlapping subcomponent densities, crucial model parameters are the 
a priori specified covariance matrix of the subcomponent means and the scale matrix of the 
inverse Wishart prior for the subcomponent covariance matrices. We select both crucial 
hyperparameters based on the variance decomposition of a mixture of mixtures model. 

We found a prior setting which is able to capture dense and connected data clusters in 
a range of benchmark data sets. However, if interest lies in detection of different cluster 
shapes, a different tuning of the prior parameters may be required. Therefore, it would be 
interesting to investigate in more detail how we can use certain prior settings to estimate 
certain kinds of data clusters. Then it would be possible to give recommendations which 
prior settings have to be used in order to capture certain types of data clusters. For instance, 
mixtures of shifted asymmetric Laplace (SAL) distributions, introduced by Franczak et al. 
(2012), have cluster distributions which are non-dense and have a strongly asymmetric 
shape with comet-like tails. In this case, the prior specifications given in Section 2 are not 
able to capture the clusters and need to be tuned to capture also this special kind of data 
clusters, see the example given in Appendix F. 

Although our approach to estimate the number of clusters worked well for many data 
sets, we encountered mixing problems with the blocked conditional Gibbs sampler outlined 
in Appendix A, in particular in high dimensional spaces with large data sets. To alleviate 
this problem, a collapsed sampler similar to Fall and Barat (2014) could be derived for 
finite mixtures. However, we leave this for future research. 
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SUPPLEMENTARY MATERIAL 


Appendix containing (A) the MCMC scheme to estimate a mixture of mixtures model, 

(B) a detailed description of the post-processing strategy based on the point process repre¬ 
sentation, (C) the simulation studies described in Section 4, (D) a description of the data 
sets studied in Section 4, (E) issues with the merging approach, and (F) estimation of data 
clusters generated by a SAL -distribution (Franczak et ah, 2012). (Appendix.pdf) 

R CODE implementing the sparse hierarchical mixture of mixtures model (Code.zip). 

A MCMC sampling scheme 

Estimation of a sparse hierarchical mixture of mixtures model is performed through MCMC 
sampling based on data augmentation and Gibbs sampling. To indicate the cluster to 
which each observation belongs, latent allocation variables S = (Si,..., Sn) taking values 
in {1,..., K} n are introduced such that 

p(yi\0i, .. .,9 K , Si = k) = Pk(y,\d k ), and Pr {Si = k\r 7 } = rj k . 

Additionally, to indicate the subcomponent to which an observation within a cluster is 
assigned to, latent allocation variables I = (/ 1; ..., In) taking values in {1,..., L} N are 
introduced such that 

Pk(yi\0k, Si = k, Ii = l) = fu{yi\pi kl , S k i) and Pr{/* = l\Si = k, w k } = w ki . 

Based on the priors specified in Section 2 . 2 , with fixed hyperparameters (eo, do, 0 ), 9o, Go, B 0 , m 0 , M 0 , u), 
the latent variables (S,I) and parameters ( 77 , w*,, fj, kl , S kh C 0 fc, b ofc , \ k j), k = 1 
I = 1,..., L, j = 1,... , r, are sampled from the posterior distribution using the following 
Gibbs sampling scheme. Note that the conditional distributions given do not indicate that 
conditioning is also on the fixed hyperparameters. 

(1) Sampling steps on the level of the cluster distribution: 

(a) Parameter simulation step conditional on the classifications S. Sample 77 |S from 
Dir(ei,e K ), e k = e 0 + N k , k = 1,, K, where N k = = k} is the 

number of observations allocated to cluster k. 
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(b) Classification step for each observation y \ conditional on cluster-specific parame¬ 
ters. For each i — 1,..., N sample the cluster assignment Si from 

Pr{Si = k\yi,d,r]} oc. r] k pk(yi\Ok), k = l,...,K, (20) 

where p k (yi\O k ) is the semi-parametric mixture approximation of the cluster den¬ 
sity: 

L 

/'/dy, ()/,) = J^«;«/v(yi|/*fcj,s w ). 

1=1 

Note that clustering of the observations is performed on the upper level of the 
model, using a collapsed Gibbs step, where the latent, within-cluster allocation 
variables I are integrated out. 

(2) Within each cluster k, k — 1,..., K: 

(a) Classification step for all observations y,, assigned to cluster k (i.e. Si = k), con¬ 
ditional on the subcomponent weights and the subcomponent-specific parameters. 
For each i e {i — 1,... r N : Si — k} sample fi from 

Pr{/i = /1yj, 6 k , Si = k} oc ww/y(y t |Mfcn / = 1, • • •, F. 

(b) Parameter simulation step conditional on the classifications I and S: 

i. Sample w fc |I, S from Dir(d k i,..., d k L ), d k i = d 0 + N kh l = 1,..., L, where 
N k i = #{/?: = l\Si = k} is the number of observations allocated to subcom¬ 
ponent l in cluster k. 

ii. For l = 1,..,, L: Sample S^/|S. I. n kl , C ok , y ~ W r (c M , C ki ), where 

c k i — Co + N k i/ 2, 

c kl = c ok + ^ fa - Pki)(yt - ■ 

i:Ii=l,Si=k 

iii. For l = 1,..., L: Sample /i kl \S, I, b ofc , S w , A fc , y ~ M r { b w , B w ), where 

Bfci = (B^ 1 + iVfc/Sfc/) -1 , 

b k i = B^^B^bofe + T, k ^N k iy k i), 
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with B ofe = \fAT k B 0 ^AT k , A k = diag(A fc i,..., X kr ), and y ki = 1/N kl y< 

i:Ii=l,Si=k 

being equal to the subcomponent mean for N k i > 0 and N k iy k i = 0, otherwise. 

(3) For each cluster k, k = 1,... ,K: Sample the random hyperparameters A k j, C 0 fc,b 0fc 
from their full conditionals: 

(a) For j = 1,..., r: Sample \ kj \h ok , n kl ,..., fx kL ~ QTQ{p kL , a kj , b kj ), where QXQ is 
the generalized inverted Gaussian distribution and 

p kL = —L/2 + v, 
a kj = 2v, 

L 

b k j 'y ^ bo k j j /Bo jj. 
i=i 

(b) Sample Cofc|Sfc 1? ..., X2 kk ~ W r (g o + Lcq, Go + Ylt=i ^h)- 

(c) Sample b ofe |B ofe , n kl ,.. ., n kL ~ jV r (m fc , M fc ), where 

Mfc = (Mq 1 + LBQ fc L ) _1 , 
m k = M fc | M 0 " L m 0 + B" 1 

\ i=i 
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B Identification through clustering in the point pro¬ 


cess representation 

Various post-processing approaches have been proposed for the MCMC output of finite or 
infinite mixture models (see, for example, Molitor et al. 2010 or Jasra et al. 2005). We 
pursue an approach which aims at determining a unique labeling of the MCMC draws 
after selecting a suitable number of clusters in order to base any posterior inference on the 
relabeled draws, such as for example the determination of cluster assignments. 

To obtain a unique labeling of the clusters, Fruhwirth-Schnatter (2006) suggested to 
post-process the MCMC output by clustering a vector-valued functional f(0 k ) of the 
cluster-specific parameters 0 k in the point process representation. The point process rep¬ 
resentation has the advantage that it allows to study the posterior distribution of cluster- 
specific parameters regardless of potential label switching, which makes it very useful for 
identification. 

If the number K of components matches the true number of clusters, it can be ex¬ 
pected that the vector-valued functionals of the posterior draws cluster around the “true” 
points ..., /(Ok)} (Fruhwirth-Schnatter, 2006, p. 96). However, in the case of an 

overfitting mixture where draws are sampled from empty components, the clustering proce¬ 
dure has to be adapted as suggested in Fruhwirth-Schnatter (2011) and described in more 
details in Malsiner-Walli et al. (2016). Subsequently, we describe how this approach can 
be applied to identify cluster-specific characteristics for the sparse hierarchical mixture of 
mixtures model. 

First, we estimate the number of non-empty clusters Kq on the upper level of the sparse 
hierarchical mixture of mixtures model. For this purpose, during MCMC sampling for each 
iteration m the number of non-empty clusters is determined, i.e. the number of clusters 

to which observations have been assigned for this particular sweep of the sampler: 

K 

i4 m) = k - H N i m) = °}, ( 21 ) 

fc=l 

where = k } i s ^ ie number of observations allocated to cluster k in 

the upper level of the mixture for iteration m and / denotes the indicator function. Then, 


following Nobile (2004) we obtain the posterior distribution of the number A' 0 of non¬ 
empty clusters Pr{A' 0 = h\y}, h = on the upper level from the MCMC output. 

An estimator of the true number of clusters A'o is then given by the value visited most 
often by the MCMC procedure, i.e. the mode of the (estimated) posterior distribution 
Pr{Ao = h\y}. 

After having estimated the number of non-empty clusters A’o, we condition the subse¬ 
quent analysis on a model with A' 0 clusters by removing all draws generated in iterations 
where the number of non-empty clusters does not correspond to A' 0 . Among the remaining 
Mq draws, only the non-empty clusters are relevant. Hence, we remove all cluster-specific 
draws 6 k for empty clusters (which have been sampled from the prior). The cluster-specific 
draws left are samples from A' 0 non-empty clusters and form the basis for clustering the 
vector-valued functionals of the draws in the point process representation into A'o groups. 

ft should be noted, that using only vector-valued functionals of the unique parameters 
6k for this clustering procedure has two advantages. First, 6k is a fairly high-dimensional 
parameter of dimension d = L — 1 + Lr(r + 3)/2, in particular if r is large, and the 
vector-valued functional allows to consider a lower dimensional problem (see also Friihwirth- 
Schnatter, 2006, 2011). In addition, we need to solve the label switching issue only on the 
upper level of the sparse hierarchical mixture of mixtures model. Thus, we choose vector¬ 
valued functionals of the cluster-specific parameters 6k that are invariant to label switching 
on the lower level of the mixture for clustering in the point process representation of the 
upper level. We found it particularly useful to consider the cluster means on the upper 
level mixture, defined by = )C[=i w ki' > 6- ki’’ ■ 

Clustering the cluster means in the point process representation results in a classification 
sequence p < - m - ) for each MCMC iteration m indicating to which class a single cluster-specific 
draw belongs. For this, any clustering algorithm could be used, e.g., A'-means (Hartigan 
and Wong, 1979) or A'-centroids cluster analysis (Leisch, 2006) where the distance between 
a point and a cluster is determined by the Mahalanobis distance, see Malsiner-Walli et al. 
(2016, Section 4.2) for more details. Only the classification sequences p which correspond 
to permutations of (1,..., AT 0 ) are used to relabel the draws. To illustrate this step, consider 
for instance, that for A' 0 = 4, for iteration m a classification sequence p^ = (1,3,4, 2) is 
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obtained through the clustering procedure. That means that the draw of the first cluster 
was assigned to class one, the draw of the second cluster was assigned to class three and so 
on. In this case, the draws of this iteration are assigned to different classes, which allows to 
relabel these draws. As already observed by Friihwirth-Schnatter (2006), all classification 
sequences p^ m \ m = 1,... , M obtained in this step are expected to be permutations, if 
the point process representation of the MCMC draws contains well-separated simulation 
clusters. 

Nevertheless, it might happen that some of the classification sequences p^ are not 
permutations. E.g., if the classification sequence p < '- rn ' t = (3,1,2,1) is obtained, then draws 
sampled from two different clusters are assigned to the same class and no unique labels 
can be assigned. If only a small fraction M 0jP of non-permutations is present, then the 
posterior draws corresponding to the non-permutation sequences are removed from the 
Mo draws with Kq non-empty clusters. For the remaining Md( 1 — Mo iP ) draws, a unique 
labeling is achieved by relabeling the clusters according to the classification sequences p^ m \ 
If the fraction M 0jP is high, this indicates that in the point process representation clusters 
are overlapping. This typically happens if the selected mixture model with Kq clusters is 
overfitting, see Friihwirth-Schnatter (2011). 

This post-processing strategy of the MCMC draws obtained using the sampling strategy 
described in Appendix A can be summarized as follows: 

1. For each iteration m — 1,..., M of the MCMC run, determine the number of non¬ 
empty clusters according to (21). 

2. Estimate the number of non-empty clusters by Kq = mode{K ^ n ^) as the value of the 
number of non-empty clusters occurring most often during MCMC sampling. 

3. Consider only the subsequence of all MCMC iterations of length M 0 where the number 

of non-empty clusters K^ n> is exactly equal to Kq. For each of the resulting m = 
1,..., Mq draws, relabel the posteriors draws Op 1 , ..., 0j^\ the weight distribution 
rii' l \ ..., as well as the upper level classifications S[ rn \ ..., S^ such that empty 
clusters, i.e. clusters with = 0, appear last. Remove the empty clusters and 

keep only the draws 6 ( \ rn \ ..., 6^ of the Kq non-empty clusters. 
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4. Arrange the Kq cluster means p ^™\..., p^f* for all Mq draws in a “data matrix” 
with Kq ■ M 0 rows and r columns such that the first Kq rows correspond to the first 
draw p^\ ..., P'p^ the next Ao rows correspond to the second draw p\\ ..., p ( ^ Q , 
and so on. The columns correspond to the different dimensions of p. Cluster all 
K 0 ■ M 0 draws into K 0 clusters using either A'-means (Hartigan and Wong, 1979) or 
A'-centroids cluster analysis (Leisch, 2006). Either of these cluster algorithms results 
in a classification index for each of the Kq ■ Mq rows of the “data matrix” constructed 
from the MCMC draws. This classification vector is rearranged in terms of a sequence 
of classifications p^ m \ m = 1,..., M 0 , where each p^ = (p[ m \ ..., P^) is a vector 
of length Kq, containing the classifications for each draw p[ m \ ..., p^' 1 at iteration 
m. Hence, p'™ 1 indicates for each single draw p^ to which cluster it belongs. 

5. For each iteration m, m — 1,..., Mo, check whether p^ is a permutation of (1,..., Kq). 
If not, remove the corresponding draws from the MCMC subsample of size Mq. The 
proportion of classification sequences of M 0 not being a permutation is denoted by 
Mo,p- 

6. For the remaining M 0 (l — M 0)P ) draws, a unique labeling is achieved by resorting the 
entire vectors of draws {0[ nL \ ..., 0 { ^} (not only p }..., P^), the weight distribu¬ 
tion Tj\ m \ ..., as well as relabeling the upper level classifications S[ m \ ..., S^ 
according to the classification sequence p( m \ 

Based on the relabeled draws cluster-specific inference is possible. For instance, a straight¬ 
forward way to cluster the data is to assign each observation y \ to the cluster S) which is 
visited most often. Alternatively, each observation yj may also be clustered based on esti¬ 
mating tik = Pr{Si = k\yi}. An estimate tik of tik can be obtained for each k — 1,..., AT, 
by averaging over Pr{S) = k\yi, given by Equation (20) using the relabeled 

draws. Each observation y* is then assigned to that cluster which exhibits the maximum 
posterior probability, i.e. 5) is defined in such a way that t i §. = maxMifc- The closer t i §. is 
to one, the higher is the segmentation power for observation i. Furthermore, the clustering 
quality of the estimated model can also be assessed based on estimating the posterior ex¬ 
pected entropy. The entropy of a finite mixture model is defined in Celeux and Soromenho 
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(1996) and also described in Friihwirth-Schnatter (2006, p. 28). Entropy values close to 
zero indicate that observations can unambiguously be assigned to one cluster, whereas large 
values indicate that observations have high a posteriori probabilities for not only one, but 
several clusters. 

To illustrate identification through clustering the draws in the point process represen¬ 
tation in the present context of a mixture of mixtures model, a sparse hierarchical mixture 
of mixtures model with K = 10 clusters and L — 4 subcomponents is fitted to the AIS data 
set (see Figure E.9 and Section 4). The point process representation of the weighted cluster 
mean draws fJL ™ ] = Yld=i °f clusters, including empty clusters, is shown in 

Figure B.4 on the left-hand side. Since a lot of draws are sampled from empty clusters, 
i.e. from the prior distribution, the plot shows a cloud of overlapping posterior distributions 
where no cluster structure can be distinguished. However, since during MCMC sampling 
in almost all iterations only two clusters were non-empty, the estimated number of clusters 
is K 0 = 2. Thus all draws generated in iterations where the number of non-empty clusters 
is different from two and all draws from empty clusters are removed. The point process 
representation of the remaining cluster-specific draws is shown in the scatter plot in the 
middle of Figure B.4. Now the draws cluster around two well-separated points, and the 
two clusters can be easily identified. 

To illustrate the subcomponent distributions which are used to approximate the cluster 
distributions, the point process representation of the subcomponent means is shown in Fig¬ 
ure B.4 on the right-hand side for the cluster discernible at the bottom right in Figure B.4 
in the middle. The plot clearly indicates that all subcomponent means are shrunken to¬ 
ward the cluster mean as the variation of the subcomponent means is about the same as 
the variation of the cluster means. 


C Simulation studies 

For both simulation studies, 10 data sets are generated and a sparse hierarchical mixture 
of mixtures model is estimated. Prior distributions and hyperparameters are specified as 
described in Section 2.1 and 2.3. MCMC sampling is run for M = 4000 iterations after a 
burn-in of 4000 draws. For the sampling, the starting classification of the observations is 
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Figure B.4: AIS data set, K = 10, L = 4, (j) B = 0.5, (j)w = 0.1: Point process representation of the 


cluster means fi k of all 10 clusters (left-hand side) and only from those where Kq — 2 (middle). Right-hand 


side: Point process representation of the means of all subcomponents forming the cluster in the bottom 
right in the plot in the middle. 

obtained by first clustering the observations into K groups using A'-means clustering and 
by then allocating the observations within each group to the L subcomponents by using 
means clustering again. The estimated number of clusters is reported in Tables C.2 
and C.3, where in parentheses the number of data sets for which this number is estimated 
is given. 

C.l Simulation setup I 

The simulation setup I consists of drawing samples with 800 observations grouped in four 
clusters. Each cluster is generated by a normal mixture with a different number of sub¬ 
components. The four clusters are generated by sampling from an eight-component normal 
mixture with component means 



1.5 6 6 1.5 8 31 31 29 


6 4 8 22.5 20 22 22 6.5 
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variance-covariance matrices 


£1 = 




S 2 


f 3.61 5.05 

V 5.05 14.44 


£4 = 


( 12.25 

V 0 




£7 = 




£s 


/ 2.25 
^4.20 



£3 


l 3.61 -5.05 

1—5.05 14.44 




1 


and weight vector 77 = l/4(l/3,1/3,1/3,1/2,1/2,1/2,1/2,1). 

In Figure C.5 the scatter plot of one data set and the 90% probability contour lines of 
the generating subcomponent distributions are shown. The first three normal distributions 
generate the triangle-shaped cluster, the next two the L-shaped cluster, and the last three 
distributions the cross-shaped and the elliptical cluster. The number of generating distri¬ 
butions for each cluster (clockwise from top left) is 1, 2, 2, and 3. This simulation setup is 
inspired by Baudry et al. (2010) who use clusters similar to the elliptical and cross-shaped 
clusters on the top of the scatter plot in Figure C.5. However, our simulation setup is 
expanded by the two clusters at the bottom which have a triangle and an L shape. Our 
aim is to recover the four clusters. 

If we estimate a sparse finite mixture model (see Malsiner-Walli et ah, 2016), which can 
be seen as a special case of the sparse hierarchical mixture of mixtures model with number 
of subcomponents L — 1, the estimated number of components is seven, as can be seen in 
the classification results shown in Figure C.5 in the middle plot. This is to be expected, 
as by specifying a standard normal mixture the number of generating normal distributions 
is estimated rather than the number of data clusters. In contrast, if a sparse hierarchical 
mixture of mixtures model with K = 10 clusters and L = 4 subcomponents is fitted to 
the data, all but four clusters become empty during MCMC sampling and the four data 
clusters are captured rather well, as can be seen in the classification plot in Figure C.5 on 
the right-hand side. 

In order to study the effect of changing the specified maximum number of clusters 
K and subcomponents L on the estimation result, a simulation study consisting of 10 
data sets with the simulation setup as explained above and varying numbers of clusters 
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L 

K 

1 

3 

4 

5 

4 

4(10) 

4(10) 

4(10) 

4(10) 

10 

7(9) 

4(10) 

4(10) 

4(10) 


6(1) 




15 

7(9) 

4(10) 

4(9) 

4(10) 


8(1) 


5(1) 



Table C.2: Simulation setup I (based on 10 data sets); true number of clusters equal to 4. Results for 
the estimated number of non-empty clusters Kq for various values of K and L. The number of data sets 
estimating the reported I\q is given in parentheses. 

K = 4,10,15 and subcomponents L = 1,3,4, 5 is performed. For each combination of K 
and L the estimated number of clusters is reported in Table C.2. 

First we study the effect of the number of specified subcomponents L on the estimated 
number of data clusters. As can be seen in Table C.2, we are able to identify the true 
number of clusters if the number of subcomponents L forming a cluster is at least three. 
I.e. by specifying an overfitting mixture with K = 10 clusters, for (almost) all data sets 
superfluous clusters become empty and using the most frequent number of non-empty 
clusters as an estimate for the true number of data clusters gives good results. If a sparse 
finite normal mixture is fitted to the data, for almost all data sets 7 normal components 
are estimated. Regarding the maximum number of clusters K in the overfitting mixture, 
the estimation results do scarcely change if this number is increased to K = 15, as can be 
seen in the last row of Table C.2. This means that also in a highly overfitting mixture, all 
superfluous clusters become empty during MCMC sampling. 

In Figure C.6, the effect of the number of subcomponents L on the resulting cluster 
distributions is studied. For the data set shown in Figure C.5, for an increasing number of 
subcomponents the estimated cluster distributions are plotted using the MAP estimates of 
the weights, means and covariance matrices of the subcomponents. The estimated cluster 
distributions look quite similar, regardless of the size of L. This robustness may be due to 
the smoothing effect of the specified hyperpriors. 


35 











Figure C.5: Simulation setup I. Scatter plot of one data set with the generating component densities 
shown with 90% probability contour lines (left-hand side), and clustering results by estimating a sparse 
hierarchical mixture of mixtures model with I\ = 10, L = 1 (middle) and K = 10, L = 4 (right-hand side). 



Figure C.6: Simulation setup I, K = 10, various values of L. For the data set in Figure C.5, the estimated 
cluster distributions (MAP estimates of means, weights, and covariance matrices of the subcomponents) 
are plotted for an increasing number of specified subcomponents L = 3,4,5 (from left to right). 

C.2 Simulation setup II 

In Section 2.3 it is suggested to specify the between-cluster variability by 0# = 0.5 and the 
between-subcomponent variability by cj)w = 0.1. As can be seen in the previous simulation 
study in Section C.l this a priori specification gives promising results if the data clusters 
are well-separated. However, in contrast to the simulation setup I, in certain applications 
data clusters might be close or even overlapping. In this case, the clustering result might be 
sensitive in regard to the specification of (j)B and 0p/. Therefore, in the following simulation 
study it is investigated how the specification of 0 b and 0yp affects the identification of data 
clusters if they are not well-separated. We want to study how robust the clustering results 
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are against misspecification of the two proportions. 

In order to mimic close data clusters, 10 data sets with 300 observations are generated 
from a three-component normal mixture, where, however, only two data clusters can be 
clearly distinguished. In Figure C.7 the scatter plot of one data set is displayed. The 300 
observations are sampled from a normal mixture with component means 


(Mi M 2 Ms) 


f 2 4.2 7.8 
1 2 4.2 7.8 


variance-covariance matrices Xi = X 2 = S 3 = l 2 and equal weights t] = (1/3,1/3,1/3). 

For various values of 0 b (between 0.1 and 0.9) and 0^/ (between 0.01 and 0.4) a sparse 
mixture of mixtures model with K = 10 clusters and L = 4 subcomponents is fitted and 
the number of clusters is estimated. For each combination of 0 b and (f>w the results are 
reported in Table C.3. 

Table C.3 indicates that if 0 b increases, also <f>w bas to increase in order to identify 
exactly two clusters. This makes sense since by increasing 0 b the a priori within-cluster 
variability becomes smaller yielding tight subcomponent densities. Tight subcomponents 
in turn require a large proportion 0^ of variability explained by the subcomponent means 
to capture the whole cluster. Thus 0w bas to be increased too. However, (ft w has to be 
selected carefully. If 0^/ is larger than actually needed, some subcomponents are likely to 
“emigrate” to neighboring clusters. This leads finally to only one cluster being estimated 
for some data sets. This is basically the case for some of the combinations of 0 b and (f>w 
displayed in the upper triangle of the table. In contrast, if (f>w is smaller than needed, 
due to the induced shrinkage of the subcomponent means toward the cluster center, the 
specified cluster mixture distribution is not able to fit the whole data cluster and two 
cluster distributions are needed to fit a single data cluster. This can be seen for some of 
the combinations of 0 b and (j>w displayed in the lower triangle of the table. 


D Description of the data sets 

The following data sets are investigated. The Yeast data set (Nakai and Kanehisa, 1991) 
aims at predicting the cellular localization sites of proteins and can be downloaded from 


37 


<t>w 

0.01 

0.1 

0.2 

0.3 

0.4 

0.1 

3(6) 

2(10) 

2(5) 

1(8) 

1(8) 


2(4) 


1(5) 

2(2) 

2(2) 

0.3 

3(6) 

2(10) 

2(8) 

2(6) 

1(7) 


2(4) 


1(2) 

1(4) 

2(3) 

0.5 

3(5) 

2(10) 

2(10) 

2(9) 

2(7) 


2(5) 



1(1) 

1(3) 

0.7 

3(7) 

2(7) 

2(10) 

2(10) 

2(10) 


2(3) 

3(3) 




0.9 

3(6) 

3(7) 

3(5) 

2(8) 

2(10) 


4(4) 

2(3) 

2(5) 

3(2) 



Table C.3: Simulation setup II (based on 10 data sets); true number of clusters equal to 2, K = 10, 
L = 4. Results for the estimated number of non-empty clusters Kq for various amounts oi (f>B and <pw- 
The number of data sets estimating the reported Kq is given in parentheses. 
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Figure C.7: Simulation setup II. Scatter plot of one data set (left-hand side), classification according to 
the generating distributions (middle) and to the clusters obtained from a mixture of mixtures with K = 10, 
L = 4, (f) B = 0.5 and (f>w = 0.1 (right-hand side). 

the UCI machine learning repository (Bache and Lichman, 2013). As in Franczak et al. 
(2012), we aim at distinguishing between the two localization sites CYT (cytosolic or 
cytoskeletal) and ME3 (membrane protein, no N-tcrminal signal) by considering a subset 
of three variables, namely McGeoch’s method for signal sequence (meg), the score of the 
ALOM membrane spanning region prediction program (aim) and the score of discriminant 
analysis of the amino acid content of vacuolar and extracellular proteins (vac). 

The Flea beetles data set (Lubischew, 1962) considers 6 physical measurements of 74 
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male flea beetles belonging to three different species. It is available in the R package 
DPpackage (Jara et al., 2011). 

The Australian Institute of Sport (AIS) data set (Cook and Weisberg, 1994) consists 
of 11 physical measurements on 202 athletes (100 female and 102 male). As in Lee and 
McLachlan (2013a), we only consider three variables, namely body mass index (BMI), lean 
body mass (LBM) and the percentage of body fat (Bfat). The data set is contained in the 
R package locfit (Loader, 2013). 

The Breast Cancer Wisconsin (Diagnostic) data set (Mangasarian et ah, 1995) describes 
characteristics of the cell nuclei present in images. The clustering aim is to distinguish 
between benign and malignant tumors. It can be downloaded from the UCI machine 
learning repository. Following Fraley and Raftery (2002) and Viroli (2010) we use a subset 
of three attributes: extreme area, extreme smoothness, and mean texture. Additionally, 
we scaled the data. 

The artificial flower data set reported by Yerebakan et al. (2014) can be downloaded 
from https://github.com/halidziya/I2GMM. It consists of 17000 two-dimensional ob¬ 
servations representing a flower shape. The data set is generated by seventeen Gaussian 
densities forming 4 clusters: nine components generate the blossom of the flower, four 
components the stem and two components each of the two leaves. Note that within each 
cluster, the generating components have the same orientation. This specification meets the 
assumption made in the infinite mixture of infinite mixtures model by Yerebakan et al. 
(2014). We used a subsample of 400 data points for our application, thus leading to the 
benchmark data sets all being of comparable size. The scatter plot of the sample is given in 
Figure D.8 on the left-hand side. If we fit a sparse mixture of mixtures model with K = 10 
clusters and L = 4 subcomponents and the usual prior settings as described in Section 2, 
the four clusters of the flower (petal, stem, and two leaves) can be clearly captured, as can 
be seen in Figure D.8, where the estimated clustering result and the corresponding cluster 
distributions are shown. 

The flow cytometry data set DLBCL contains intensities of markers stained on a sample 
of over 8000 cells derived from the lymph nodes of patients diagnosed with Diffuse Large 
B-cell Lymphoma (DLBCL). The aim of the clustering is to group the individual cell data 
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Figure D.8: Flower data set. Boxplot of a sample with 400 data points (left-hand side), the estimated 
clusters for K = 10, L = 4, (/>b = 0.5, <f>w = 0.1, = 10 (middle), and the corresponding cluster 

distributions (right-hand side). 


measurements into only a few groups on the basis of similarities in light scattering and 
fluorescence, see Aghaeepour et ah (2013) for more details. For this data set, class labels of 
the observations partitioning the data into four classes are available which were obtained by 
manual partitioning (“gating”). The data set is available in the R package EMMIXuskew 
(Lee and McLachlan, 2013b) as data set DLBCL with the corresponding class labels in 
true.clusters. 

The flow cytometry data set GvHDBOlcase by Brinkman et al. (2007) consists of 12442 
six-dimensional observations which represent a blood sample from a subject who developed 
Graft versus Host disease (GvHD). GvHD is a severe complication following a blood and 
marrow transplantation, when donor immune cells in the graft attack the body cells of the 
recipient. The data were analyzed first by Brinkman et al. (2007). Lo et ah (2008) fitted 
a Student-f mixture model to this data and estimated 12 clusters using the EM algorithm. 
In the Bayesian framework, Friihwirth-Schnatter and Pyne (2010) fitted finite mixtures 
of skew-normal and skew-f distributions and found 12 and 9 clusters. By comparing this 
sample to a control sample from a patient who had a similar transplantation but did not 
develop the disease, Brinkman et al. (2007) found a very small cluster of live cells (high 
“FSC”, high “SSC”) in the sample with a high expression in the four markers (“CD4+”, 
“CD8/3+”’, “CD3+”, “CD8+”). This cluster was not present in the control sample and 
seems to be correlated with the development of GvHD. 
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For the data sets with known class labels, the clustering result of the estimated models 
is measured by the misclassihcation rate and the adjusted Rand index (Hubert and Arabie, 
1985). To calculate the misclassihcation rate of the estimated model, the “optimal” match¬ 
ing between the estimated cluster labels and the true known class labels is determined as 
the one minimizing the misclassihcation rate over all possible matches for each of the sce¬ 
narios. The misclassihcation rate is measured by the number of misclassihed observations 
divided by all observations and should be as small as possible. 

The adjusted Rand index (Hubert and Arabie, 1985) is used to assess the similarity 
between the true and the estimated partition of the data. It is a corrected form of the 
Rand index (Rand, 1971) which is adjusted for chance agreement. An adjusted Rand index 
of 1 corresponds to perfect agreement of two partitions whereas an adjusted Rand index 
of 0 corresponds to results no better than would be expected by randomly drawing two 
partitions, each with a hxed number of clusters and a hxed number of elements in each 
cluster. 


E Issues with the merging approach 

The merging approach, which consists of hrst fitting a finite mixture of Gaussians to suit¬ 
ably approximate the data distribution and subsequently combines components to clusters, 
is susceptible to yield poor classihcations, since the resulting clusters can only emerge as 
the union of components that have been identified in the previous step. For illustration, the 
AIS data (see Appendix D) are clustered using function clustCombi (Baudry et ah, 2010) 
from the R package mclust (Fraley et al., 2012). The results are shown in Figure E.9. The 
first step results in a standard Gaussian mixture with three components (left-hand plot), 
and subsequently all data in the smallest component are merged with one of the bigger 
components to form two clusters (middle plot) which are not satisfactorily separated from 
each other due to the misspecihcation of the standard Gaussian mixture in the hrst step. 
In contrast, the sparse hierarchical mixture of mixtures approach we develop in the present 
paper identifies two well-separated clusters on the upper level of the hierarchy (right-hand 
plot). 
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Figure E.9: AIS data set, variables “X.Bfat” and “LBM”. Scatter plots of the observations with different 
estimated classifications based on Mclust (left-hand side), combiClust (middle), and the sparse hierarchical 
mixture of mixtures approach developed in this paper ( I\ = 10, L = 4) (right-hand side). 

F Fitting a mixture of two SAL distributions 

Although it is not the purpose of our approach to capture non-dense data clusters, we 
apply it to the challenging cluster shapes generated by shifted asymmetric Laplace (SAL) 
distributions, which are introduced by Franczak et al. (2012) in order to capture asymmetric 
data clusters with outliers. We sampled data from a mixture of two SAL distributions 
according to Section 4.2 in Franczak et al. (2012). The data set is shown in Figure F.10 
on the left-hand side. 

If we fit a sparse hierarchical mixture of mixtures model with K = 10 clusters and L = 4 
subcomponents and priors and hyperpriors specified as in Sections 2.1 and 2.3, four clusters 
are estimated, as can be seen in the middle plot of Figure F.10. Evidently, the standard 
prior setting, tuned to capture dense homogeneous data clusters, performs badly for this 
kind of clusters. Thus, in order to take the specific data cluster shapes into account, we 
adjust the prior specifications accordingly. A data cluster generated by a SAL distribution 
is not homogeneously dense, it rather consists of a relatively dense kernel on one side of the 
cluster and a non-dense, light and comet-like tail with possibly extreme observations on 
the other side. Therefore within a cluster, subcomponents with very different covariance 
matrices are required in order to fit the whole cluster distribution. Since specification of 
hyperpriors on A kj and Cofc has a smoothing and balancing effect on the subcomponent 
densities, we omit these hyperprior specifications, and choose fixed values for k = 1 ,..., K , 
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Figure F.10: Samples from a mixture of two SAL distributions (left-hand side), the estimated clusters 
for K = 10 ,L = 4, cj)B = 0.5, fw = 0.1, = v 2 = 10 (middle), and for K = 10, L = 5, (j>s = 0.4, fw = 0.2, 

with fixed hyperparameters C 0 k = go • G^ 1 and A k j = 1 (right-hand side). 

i.e. C 0 fc = go ■ G 0 1 and X kj = 1 , j = 1 ,. .., r. 

Additionally, in order to reach also extreme points, we increase both the number of 
subcomponents to L = 5 and the a priori variability explained by the subcomponent 
means to (f>w = 0.2. At the same time we adjust the proportion of heterogeneity explained 
by the cluster means by decreasing 0 b to 0.4, thus keeping the subcomponent covariance 
matrices large. If we estimate again a sparse hierarchical mixture of mixtures model with 
these modified prior settings, the two clusters can be identified, see Figure F.10 on the 
right-hand side. 
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